function remove_staves(sheet,staffline_height, staffspace_height)

DEBUG = 1;

% calculate distance transform of original image
distance_transform = bwdist(sheet,'euclidean'); % distance transform of original image
% calculate image skeleton
skeleton = bwmorph(~sheet,'skel',Inf);

% remove spurs that are shorter than staffline_height/2
for i = 1:floor((staffline_height+1)/2)
    skeleton = bwmorph(skeleton, 'spur');
end

% split skeleton
allsegs = split_skeleton(skeleton, distance_transform, ...
    staffline_height, staffspace_height);

% filter staff segment candidates
[candidates, nocandidates] = filter_stafsegments(allsegs, ...
    distance_transform, staffline_height, staffspace_height);

if DEBUG
    [r, c] = size(skeleton);
    rgb = zeros(r, c, 3);
    rprops = regionprops(allsegs,'PixelIdxList');
    temp_image = zeros(r, c);
    for i = 1:length(candidates)
        temp_image(rprops(candidates(i)).PixelIdxList) = 255;
    end
    rgb(:,:,2) = temp_image(:,:);
    
    temp_image = zeros(r, c);
    for i = 1:length(nocandidates)
        temp_image(rprops(nocandidates(i)).PixelIdxList) = 255;
    end
    rgb(:,:,1) = temp_image(:,:);
    
    imwrite(rgb,'debug.png','png');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function allsegs = split_skeleton(skeleton, distance_transform, staffline_height, staffspace_height)
% find branch points
branch_points=bwmorph(skeleton, 'branchpoints');

% find corner points with an angle below 135 degrees
x = im2single(skeleton); 
hcornerdet = video.CornerDetector( ... 
   'Method', ...,
   'Local intensity comparison (Rosen & Drummond)', ... 
   'CornerLocationOutputPort', false, ... 
   'MetricMatrixOutputPort', true, ... 
   'MaximumAngleThreshold', '135.0');  
 corners = step(hcornerdet, x);  
 corners = ~(~corners); % convert to logical

% Create shape inserter object
circle = video.ShapeInserter('Shape', 'Circles', ...
    'Fill', true, ...
    'FillColor', 'White', ...
    'Opacity', 1);

% Create coordinates of shapes for branch points
[i, j] = find(branch_points);
pts = zeros(3, length(i));
for k = 1:length(i)
    pts(:,k) = [i(k), j(k), distance_transform(i(k), j(k)) ];
end

% insert shapes into branch points
branch_points = step(circle, double(branch_points), pts);
branch_points = ~(~branch_points); % convert to binary

% Create coordinates of shapes for corner points
[i, j] = find(corners);
pts = [i, j, max([floor((staffspace_height+1)/2), 5]) .* ones(length(i), 1)];
pts=transpose(pts);

% insert shapes into corner points
corners = step(circle, double(corners), pts);
corners = ~(~corners); % convert to binary

% unify branch and corner points
splitting_points = or(corners, branch_points); % splitting points that will be removed

% Remove splitting points
skeleton = and(skeleton, not(splitting_points));

% splitted skeleton is connected component labeled
allsegs = bwconncomp(skeleton);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [candidates, nocandidates] = filter_stafsegments(allsegs, distance_transform, staffline_height, staffspace_height)

rprops = regionprops(allsegs,'BoundingBox','PixelList','Centroid','Area');
candidates = [];
nocandidates = [];
maxstraightness = max([0.125*staffline_height^2, 3]);
distmax = staffline_height;
distmin = max([staffline_height/4,1]);

for i = 1:allsegs.NumObjects
    if rprops(i).BoundingBox(3) > rprops(i).BoundingBox(4) % if bounding box of the segment wider than tall
        y=rprops(i).PixelList(:,2); % Vertical coordinates of the segment
        x=rprops(i).PixelList(:,1); % Horizontal coordinates of the segment
        [b,bint,r,rint,stats] = regress(y,[ones(size(x)) x]);
        orientation_angle = 180*atan(b(2))/pi; % calculate orientation angle of the least square fitted line of the segment in degrees
        straightness = stats(4); % mean square deviation from least square fitted line
        if (abs(orientation_angle) < 25) && ...
                (straightness < maxstraightness ) && ...
                distance_percentage(rprops(i).PixelList, ...
                distance_transform, ...
                distmin, distmax) > 0.7
            candidates = [candidates i];
        end
    else
        nocandidates = [nocandidates i];
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function percentage = distance_percentage(points ,distance_transform, distmin, distmax)

[r, c] = size(points);
counter = 0;
for i = 1:r
    distance = distance_transform(points(i,2), points(i,1));
    if distance < distmax && distance > distmin
        counter = counter + 1;
    end
end

percentage = counter / r;
            